/* File converted with FORTRAN CONVERTER utility.
   FORTRAN CONVERTER is written by Grigoriev D., 2081/4.*/

/*MARKERS: //! = probably incorrect, //? = possibly incorrect*/

#include "forsythe.h"


SPLINE :: SPLINE(int _N, REAL *_X, REAL *_Y)
{

//      ВЫЧИСЛЯЮТСЯ КОЭФФИЦИЕНТЫ В(Г), C(I) И D(I), I=1,2,..., N, ДЛЯ КУБИЧЕСКОГО ИНТЕРПОЛЯЦИОННОГО СПЛАЙНА
//      S(X) = Y(I)+ B(I)*(X-X(I))+C(I)*(X-X(I))**2+H-D(I)*(X-X(I))**3
//      ДЛЯ X(I) .LE. X .LE. X(H-l)
//
//      ВХОДНАЯ ИНФОРМАЦИЯ..
//      N = ЧИСЛО ЗАДАННЫХ ТОЧЕК ИЛИ УЗЛОВ (N .GE. 2)
//      Х = АБСЦИССЫ УЗЛОВ В СТРОГО ВОЗРАСТАЮЩЕМ ПОРЯДКЕ
//      Y = ОРДИНАТЫ УЗЛОВ
//
//      ВЫХОДНАЯ ИНФОРМАЦИЯ...
//      В, С, D = МАССИВЫ ОПРЕДЕЛЕННЫХ ВЫШЕ КОЭФФИЦИЕНТОВ СПЛАЙНА.
//      ЕСЛИ ОБОЗНАЧИТЬ ЧЕРЕЗ Р СИМВОЛ ДИФФЕРЕНЦИРОВАНИЯ, ТО
//      Y(I)=S(X(I))
//      B(i)=sP(X(i))
//      C(I)= SPP(X(l))/2
//      D(I) = SPPP(X(I))/6 (ПРАВОСТОРОННЯЯ ПРОИЗВОДНАЯ)
//      С ПОМОЩЬЮ СОПРОВОЖДАЮЩЕЙ ПОДПРОГРАММЫ-ФУНКЦИИ SEVAL МОЖНО ВЫЧИСЛЯТЬ ЗНАЧЕНИЯ СПЛАЙНА.

N = _N;

X = new REAL[N];
Y = new REAL[N];
B = new REAL[N];
C = new REAL[N];
D = new REAL[N];

for(int i = 0; i < N; i++)
{
    X[i] = _X[i];
    Y[i] = _Y[i];
}

int NM1,IB,I;
REAL T;

NM1=N-1;
if(N<2)return;  //?
if(N<3)goto _50;

//     ПОСТРОИТЬ ТРЕХДИАГОНАЛЬНУЮ СИСТЕМУ ПОДПРОГРАММЫ SPLINE И SEVAL
//     В =ДИАГОНАЛЬ, D = НАДДИАГОНАЛЬ, С =ПРАВЫЕ ЧАСТИ

D[0]=X[1]-X[0];
C[1]=(Y[1]-Y[0])/D[0];
for( I=2; I<=NM1; I++){   //? target=10
D[I-1]=X[I]-X[I-1];
B[I-1]=2.*(D[I-2]+D[I-1]);
C[I]=(Y[I]-Y[I-1])/D[I-1];
C[I-1]=C[I]-C[I-1];
//_10:;
}                         	// CONTINUE
//     ГРАНИЧНЫЕ УСЛОВИЯ. ТРЕТЬИ ПРОИЗВОДНЫЕ В ТОЧКАХ
//     X (I) И X (N) ВЫЧИСЛЯЮТСЯ с ПОМОЩЬЮ РАЗДЕЛЕННЫХ РАЗНОСТЕЙ

B[0]=-D[0];
B[N-1]=-D[N-2];
C[0]=0.;
C[N-1]=0.;
if(N==3)goto _15;
C[0]=C[2]/(X[3]-X[1])-C[1]/(X[2]-X[0]);
C[N-1]=C[N-2]/(X[N-1]-X[N-3])-C[N-3]/(X[N-2]-X[N-4]);
C[0]=C[0]*D[0]*D[0]/(X[3]-X[0]);
C[N-1]=-C[N-1]*D[N-2]*D[N-2]/(X[N-1]-X[N-4]);

//     ПРЯМОЙ ХОД

_15:;
for( I=2; I<=N; I++){   //? target=20
T=D[I-2]/B[I-2];
B[I-1]=B[I-1]-T*D[I-2];
C[I-1]=C[I-1]-T*C[I-2];
//_20:;
}                         	// CONTINUE

//     ОБРАТАНЯ ПОДСТАНОВКА

C[N-1]=C[N-1]/B[N-1];
for( IB=1; IB<=NM1; IB++){   //? target=30
I=N-IB;
C[I-1]=(C[I-1]-D[I-1]*C[I])/B[I-1];
//_30:;
}                         	// CONTINUE

//     В C(I) ТЕПЕРЬ ХРАНИТСЯ ВЕЛИЧИНА SIGMA(I), ОПРЕДЕЛЕНC НАЯ В § 4.4.
//
//    С ВЫЧИСЛИТЬ КОЭФФИЦИЕНТЫ ПОЛИНОМОВ

B[N-1]=(Y[N-1]-Y[NM1-1])/D[NM1-1]+D[NM1-1]*(C[NM1-1]+2.*C[N-1]);
for( I=1; I<=NM1; I++){   //? target=40
B[I-1]=(Y[I]-Y[I-1])/D[I-1]-D[I-1]*(C[I]+2.*C[I-1]);
D[I-1]=(C[I]-C[I-1])/D[I-1];
C[I-1]=3.*C[I-1];
//_40:;
}                         	// CONTINUE
C[N-1]=3.*C[N-1];
D[N-1]=D[N-2];
return;  //?

_50:;
B[0]=(Y[1]-Y[0])/(X[1]-X[0]);
C[0]=0.;
D[0]=0.;
B[1]=B[0];
C[1]=0.;
D[1]=0.;
return;  //?
}                              	// END

SPLINE :: ~SPLINE()
{
    delete [] X;
    delete [] Y;
    delete [] B;
    delete [] C;
    delete [] D;
}

REAL SPLINE :: Eval(REAL U)
{


//  ЭТА ПОДПРОГРАММА ВЫЧИСЛЯЕТ ЗНАЧЕНИЕ КУБИЧЕСКОГО СПЛАЙНА
//
//  SEVAL = Y(I)4-B(I)*(U—X(I))4-C(I)*(U—Х(1))**2+ D(I)*(U—Х(1))**3
//  ГДЕ X(I) .LT. U .LT. Х(14-1). ИСПОЛЬЗУЕТСЯ СХЕМА ГОРНЕРА
//  ЕСЛИ U .LT. Х(1), ТО БЕРЕТСЯ ЗНАЧЕНИЕ 1 = 1.
//  ЕСЛИ U .GE. X(N), ТО БЕРЕТСЯ ЗНАЧЕНИЕ I = N,
//
//  ВХОДНАЯ ИНФОРМАЦИЯ..
//  N = ЧИСЛО ЗАДАННЫХ ТОЧЕК
//  U = АБСЦИССА, ДЛЯ КОТОРОЙ ВЫЧИСЛЯЕТСЯ ЗНАЧЕНИЕ СПЛАЙНА
//  X, Y = МАССИВЫ ЗАДАННЫХ АБСЦИСС И ОРДИНАТ
//  В, С, D =МАССИВЫ КОЭФФИЦИЕНТОВ СПЛАЙНА, ВЫЧИССЛЕННЫЕ ПОДПРОГРАММОЙ SPLINE
//  ЕСЛИ ПО СРАВНЕНИЮ С ПРЕДЫДУЩИМ ВЫЗОВОМ U НЕ
//  НАХОДИТСЯ В ТОМ ЖЕ ИНТЕРВАЛЕ, ТО ДЛЯ РАЗЫСКАНИЯ
//  НУЖНОГО ИНТЕРВАЛА ПРИМЕНЯЕТСЯ ДВОИЧНЫЙ ПОИСК.

int J,K;
REAL DX;
static int I=1;
if(I>=N)I=1;
if(U<X[I-1])goto _10;
if(U<=X[I])goto _30;

//  ДВОИЧНЫЙ ПОИСК

_10:;
I=1;
J=N+1;
_20:;
K=(I+J)/2;
if(U<X[K-1])J=K;
if(U>=X[K-1])I=K;
if(J>I+1)goto _20;

//  ВЫЧИСЛИТЬ СПЛАЙН

_30:;
DX=U-X[I-1];
return Y[I-1]+DX*(B[I-1]+DX*(C[I-1]+DX*D[I-1]));
}                         	     	// END
